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SMM-UVSP OZONE PROFILE INVERSION PROGRAMS 


1.0 INTRODUCTION 

The documentation for the computer programs used to process the Ultra- 
Violet Spectrometer Polarimeter (UVSP) data taken by the Solar Maximum Mission 
satellite (SMM) to obtain the profiles of ozone in solar occultation is pre- 
sented In this report. The next section presents an outline of the experi- 
ments performed together with a brief description of the technique used to 
invert the occultation data. Section 3 describes the main features of the 
programs and gives a pictorial *iew of how they link together. The final sec- 
tion presents a brief discussion of the results of this effort. Each of the 
individual programs and subprograms is described in Appendix 1 with the flow 
chart and a listing. Appendix 2 provides a user's manual with complete in- 
structions for the use of the programs. Appendix 3 contains a description of 
the geometry of the occultation and how it is computed. 


2.0 DESCRIPTION OF EXPERIMENT 

The UltraVi-olet Spectrometer Photometer (UVSP) aboard the Solar 
Maximum Mission (SMM) satellite was used to obtain profiles of the atmospheric 
ozone in occultation. The instrument was programmed to use a given slit and 
to measure the magnitude of the solar flux that was transmitted through the 
atmosphere as the sun was appearing or disappearing behind the earth. A more 
complete description of the experiment is given in Reference 1. 

In order to perform the inversion of the occultation data to obtain 
the ozone profile, two pieces of data must be used. The first is obviously 
the occultation profile itself; the second is the satellite ephemeris. The 
SMM was not capable of including a description of its location on the down- 
linked data stream sufficiently accurate for the purposes of inversion. There- 
fore, it was necessary to obtain the ephemeris from other sources namely from 
the satellite emphemeris group at Goddard Space Flight Center. The ephemeris 
must then be melded with the occultation profile with care to ensure that the 
two independent time streams match correctly.. This in fact turned out to be 
the most difficult part of the effort. 

The basic method used for the inversion of the occultation data 
starts from the relation 


x *ex 
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where is the transmittance through the tangent path, is the ozone cross 
section at the wavelength of operation , and [Oj] is the ozone concentration 
profile. The integral is to be taken over half the tangent path, that is from 
the tangent height to infinity or in practice, to some high altitude to be 
determined. This method of inversion has been discussed several times and we 
refer the interested reader to References 1 and 2. Equation 1 may be recast 
to a more tractable format 

/[ 0 3 ]dl-ln(l/T)/2o x (2) 

This is a relatively simple linear integral equation and may be solved by stan- 
dard matrix methods. It may also be pointed out that this method is entirely 
equivalent to a linear multiple regression for the concentration profile. 

3.0 OVERALL DESCRIPTION OF THE INVERSION PROGRAMS 

The inversion algorithms used for this effort were developed and run 
on a Digital Equipment Corporation (DEC) computer, a PDP-11/23 running the 
UNIX operating system. The programs were all written in Fortran 77 but were 
written such that a-~Fortran IV compiler with file commands (OPEN and CLOSE) 
could also be used. In fact the initial development was done on the DEC 
PDP-11/03 running RT-11 and using the DEC Fortran IV compiler. For those with- 
out a Fortran 77 compiler we describe in the appropriate places vrfiat few and 
minor changes must be made to run these programs. 

The SMM UVSP occultation data used in this effort was provided on 
nine track tape in standard DEC format. The files were binary data with 512 
bytes or 256 PDP-11 words per block. The first block of a file for a given 
experiment contained information about the experiment such as the date, time, 
wavelength, etc. The second block was a header block for the next logical re- 
cord which consisted of sixteen file blocks, containing the digitized output 
of the counters measuring the solar flux as observed by the system. If there 
was more than one logical record for the experiment of interest, each such log- 
ical record of sixteen blocks of counter data was preceded by a header block 

as well. The header blocks contained such information as the start time of 

the record, end time, etc. A more complete description of the format of the 

UVSP final data tapes will be given in the description of the programs that 

read them. 



The ephemeris data was also provided on nine track tape in DEC 
format. This data consisted of double precision floating point numbers for 
the most part with the remainder being integer - descriptions of dates and 
times, etc. In order to use the ephemeris data in an optimun manner for the 
UVSP inversion it was necessary to decrease somewhat the time precision of the 
data. The data is given at one second intervals. The UVSP data is presented 
with a much finer time resolutions and it is necessary to provide some form of 
interpolation on the ephemeris in order to obtain the satellite position accur- 
ately at an arbitrary time. It was found that the emphemeris time intervals 
of one second was too fine for convenience since a third order spline fit to 
the ephemeris was selected as the most appropriate interpolation method. 

In Figure 1, we give an overview of the set of programs used for the 
inversion. They may be easily be combined by using any of the usual command 
string interpreters on modem operating systems such as the UNIX Shell or 
RT-ll's CSI. The following paragraphs describe the basic features of these 
programs and their interactions with each other. The detailed descriptions 

are left for Appendix 2. 

The general _flow of the analysis is shown in Figure 1. The two data 
streams are shown at the top of the Figure as files with extension *.AQ for 
the ephemeris data and *.FD for the occultation data. (Note that the * is the 
usual DEC convention for a "wildcard" filename, i.e. it stands for any name.) 
These are processed by the data handling programs AQTOPT and DATPGM respect- 
ively in order to produce more convenient data files concentrated on the time 
frame of the occultation itself. These files are given the files extensions 

*.PT and *.0C respectively. The program TANSFM then reads both files to pro- 

duce another data file which contains the occultation profiles as a function 
of the acutual tangent point location expressed as a vector with components 
height, latitude and longitude. This is written on a file with extension 

*.TH. The latter file contains more data than needed for the analysis of Equa- 
tion 2. Thus the program BINS is used to compress and smooth the data Into a 
more convenient format on the files *.BN. Finally, the program REAL03 per- 
forms the required inversion based on Equation 2. 

The program AOTOPT was written to read the ephemeris tape and write a 
new file containing only those times around the occultation period and at ten 
second intervals. In order to be sure that the spline Interpolation worked 
properly, its output at ten seconds was compared with the same for five second 
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intervals for a number of selected cases. No appreciable difference was detec- 
ted. Thus we concluded that the satellite ephemeris is well described by a 
spline fit at ten second Intervals. 

In practice, It is more convenient to read all the ephemeris files 
off the tapes and store them on disk. We have done so and give the files the 
DEC file extension *.A0 where * Is the usual DEC wildcard convention for the 
file name. Note that each of the UVSP experiments was given a sequential 
number with the file name starting with the letter "V", e.g. "V00513". Thus 
the ephemeris data for this experiment would have the file name "V00513.A0". 
The program "AOTOPT" creates the shorter ephemeris file "*.PT" following the 
convention just described. 

The UVSP "final data" files (*.FD) were also read from tape to disk. 
The program DAT PGM provides a listing of the experiment and block header in 
printed format. It also compresses the file to contain only the data In two 
blocks around the occultation on the file *.0C. The user is asked for the 
block location containing the start of the occultation. This must be known 
before the program is run. Another program (DATFND) based on the structure of 
DATPGM, was used to enquire for these quantities by looking through the *.FD 
files. It would be possible to automate this procedure but this was not found 
necessary for the present amount of data. The data for the occultation frames 
is then written on the files *.0C for use by the progrms TANSMM and BINS. 

The program TANSMM reads the ephemeris files *.PT and performs a 
cubic spline fit to each of the three Cartesian coordinates of the orbit. The 
spline is then used to interpolate to get a precise position of the sensor at 
each of the times at which a count was recorded. The timing data Is read by 
TANSMM for the files *.0C as written by DATPGM. The solar position is obtain- 
ed by using a Chebychev polynomial fit to the solar ephemeris for the year 
1980 as published by the U.S. Naval Observatory^. This figure of the earth 

is also approximated by an ellipsoid of revolution 1 - . It was found that this 
simple form was entirely adequate for current purposes since the excursions of 
the measured figure of the earth from this figure were smaller than the preci- 
sion of the experiment. The geometry of the line of sight was determined from 
the pitch, roll and yaw of the satellite and from the location of the center 
of the solar disk. A series of rotation matrices was used to transform from 
one coordinate system to another. A simple vector relation was then used to 
determine the tangent height location. 
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The program BINS was used to perform a seven point running average to 
smooth the data. In addition, the amount of data was cut to Include only the 
region immediately around the occultatlon. For convenience In using the rela- 
tively small address space of the PDP-11, the data was further compressed by 
performing a binning of the data. The number of points to be averaged Into 
each bln Is a user input option. In practice it was found that three was a 
good number to use for this purpose. The output of these procedures was writ- 
ten on files *.BN. 

The final program used in this procedure is called REAL03. It reads 
the binned data files *.BN and performs the inversion following the algorithm 
for the solution of Equation 2. The output is printed and no output disk file 
Is made. For this purpose, we used the UNIX redirection features to good 
effect . 

4.0 DISCUSSION 

The programming effort described in this report shows the power of 
modern small computers which may be used for purposes for which only a large 
mainframe' could suffice only a few years ago. Indeed most of the work describ- 
ed here could have been carried out on many of the latest generation of person- 
al computers. The only problem would have been the use of the nine track tape 

drive for the data files. Even this could have been resolved easily by using 

a machine with a tape drive and down- loading the required data over a communi- 

cations link. 

The inversion technique employed allows a very simple algorithm to be 
used. A more complicated situation such as several species absorbing at the 
wavelength of Interest would have been quiet a bit more complicated but still 
tractable through the use of least squares methods. The main limitation 
arises from the relatively small address space allowed by the PDP-11. It 
might be pointed out that some of the new 16-bit microprocessors have a much 
larger address space and this limitation would not apply. 
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APPENDIX A 


PROGRAM DESCRIPITONS , 

In this appendix we present the write-up and descriptions of each of 
the programs and subprograms used In the ozone Inversion. An attempt has been 
made to follow the order of use In performing an Inversion of a given set of 
occultatlon data. 

1. PROGRAM DATFND 

This program Is used to search through the blocks of the given final 

data file Vxxxxx.FD., where Vxxxxx Is the experiment number of the occultatlon 

of interest. (Each UVSP experiment had been assigned such a number as It was 
accepted for Inclusion In the dally SW experiment program). The user is 
prompted to enter this number namely "Vxxxxx"; for example It might be V01815. 
The number Is then concatenated with the file extension ".FD" and the corres- 
ponding final data file is opened. The user Is prompted again for the number 
of a block to be examined as to whether it contains the beginning of the occul- 
tation or not. The answer to this question is immediately apparent when the 
output is listed on the crt screen or printed on a hardcopy terminal. If it 

is the proper block the user is to make a note of it and to either end the ex- 

ecution of DATFND or to request another experiment file. If it is not an oc- 
cultation block, one needs to input another candidate block number and contin- 
ue the search until the occultation block is found. Note that a UVSP block 
consists of 256 sixteen bit integers and an occultation can extend over two 
such blxks but not more. Thus what one wants is the nubmer of the first of 
these two blocks. 

As pointed out in the main body of the text it would be possible to 
program this procedure but it was not felt worthwhile for the current amount 
of SMM data. When the .9<M repaired in-orbit, the amount of data should un- 
dergo a dramatic increase. It will be necessary then to provide an automatic 
procedure for detecting the start of the occultation. 

2. PROGRAM DAT PGM 


The program DATPGM is used to read the file header block and the 


"page" header block to obtain the pointing and timing Information for the 
given experiment. The occultatlon block number Is then used to locate the two 
contiguous blocks which contain the actual occultat16n profile. 

The user Is prompted for the file specification which Is the experi- 
ment number Vxxxxx discussed In the previous program description for DATFND. 
The experiment number read Is then concatenated with the final data extension 
".FD." The next prompt Is for the "header block number" to which one replies 
with an Integer giving the block number of the "page" (Ihblck). One Is then 
prompted to give the number of the actual block (Iblock) on which the occulta- 
tlon starts (found by using DATFND) and the nubmer of the actual page ( 1 page) . 

It should be noted that there Is one page header block for every sixteen 
blocks of UVSP count data. Thus the very first header block would have 
1hb1ck>2 and lpage-1. The second header block would have 1 hbl ck- 18 and 
1page-2 and so forth. The .FD file Is then opened and the very first record 
(the file header record) Is read. Various quantities are printed out for the 
user's reference and some floating point numbers are computed, for example the 
experiment .start and stop times and the wavelength defining the spectrometer 
scan. The page header block is then read and more quantities are printed for 
reference and the pointing angles (pitch, roll and yaw) are computed together 
with start and stop times of the page. Note that these quantities have been 
read into two different arrays Ismm and ismn4. The first array Is for 2 byte 
Integers and the second Is for 4 byte Inters. 

The actual blocks containing the occultatlon are then read Into 
arrays 1smm(256) and ismml(256) and, b/ an equivalence statement, the full 
array of the two blocks containing the occultation profile, 1smmfl(512), Is 
obtained simultaneously. The final data file Is then closed. 

The experiment number is concatenated with the string ".OC" and a 
file with this name Is created as a sequential formatted file. Note that » ,ii s 
file may be printed for the user's convenience If he so chooses. The floating 
point array data(512) Is created by setting It equal to the array 1smmfl(512) 
in a do loop. The maximum value of the array data(512) Is then found by a 
call to FUNCTION RMXMN and the array data(512) Is normalized using this maxi- 
mum value. The values of the timing and angular Information needed for the 
Inversion are then written to the .X file followed by the array data. 
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3. PROGRAM AOTOPT 

I This program reads the ephemerls tape and prepares the shortened form 

of the ephemerls for use In the spline fit. The File names are obtained In 

the same way as outlined for DATFND and OATPGM. The ephemerls has the exten- 
sion H .PT" and the short ephemerls has the extension M .PT". The user Is prcmp- 

* ted for the start and stop times of Interest In seconds of day. The .AO file 

Is then read and for every tenth second the Cartesian coordinates are written 

to the .PT file together with the actual time In seconds. For convenience 

some output Is printed as the procedure Is followed. This may be kept as 

* record of the process. 

4. PROGRAM TANSMM 

* This Is the most complicated of the programs discussed In this re- 
port. Its purpose Is to read the two files *.PT and *.0C to obtain the pro- 
file of the data, not as a function of time but as a function of the tangent 
height of the ray path to the point on the sun that was being observed through 

>' the spectrometer slit. In the given experlement. As mentioned In the main 

text, to do this It Is necessary to have an accurately matched set of the two 
time series, namely the UVSP occultatlon data and the satellite ephemerls. 

The user first supplies the familiar experiment number and this Is 
concatenated with the .OC extension. The timing and anglular Information Is 
then read off this file and It is closed. The .PT file Is then opened and the 
ephemerls data Is read by calling subroutine EPHTST. This routine performs 
the calculation of the spline coefficients by calling the SPLINE subroutine. 
More details will be given in the appropriate section. 

The pointing Information Is given with respect to the solar axis of 
rotation and thus this direction must be calculated and stored In the array 
solaxs(3). Seven of the .FD files we received had timing errors on them and 
an offset of 262.144 secs had to be subtracted. This Is done In a do loop 
where a test is made to ascertain whether the present case Is one of the seven. 

The main loop of the program over each of the 512 data points contain- 
ed In the occultatlon Is then entered. The first action taken Is to compute 

the Cartesian position of the sun center by a call to the solar ephemerls sub- 

routine SOLEPH. The Cartesian satellite position at the current time Is then 
computed by using the cubic spline evaluating function SEVAL. The radial posl- 

<■ 14 J 
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tion of the satellite is then computed by a call to the subroutine LATLON, and 
the unit vector from the satellite to the solar center is obtained by a call 
to the subroutine UNIVEC. The rotation matrices needed are then calculated in 
subroutine ROLLIT, and the tangent height of the satellite- solar center path 
is obtained from a call to subroutine TANHGT. The first rotation is performed 
by a call to subroutine TANHGT. The first rotation is performed by a call to 
subroutine ROTMTX and this is followed iy a translation. The next rotation 

precedes the calculation of the tangent point coordinates. The transformation 
process is then reversed by two calls to subroutine ROTMTX followed by a trans- 
lation in the reverse sense to the preceding one and the final rotation is 
then processed. The calculation of the actual tangent point coordinates for 
the current sight path in subroutine TANHGT is preceded by the determination 
of the appropriate unit vector in subroutine UNIVEC. Finally, the figure of 
the earth is used to compute an accurate value for the height of the tangent 
point. To provide some kind of record of the computation the values for every 
tenth point obtained are printed. 

When the loop has finished a new file is created with the extension 
“.TH" as a sequential, formatted file which may be printed if desired. The 
quantities printed include the latitude and longitude of the last tangent 
point computed and the array of tangent heights, htan(512). This file is then 
closed and the program ends. 

5. PROGRAM BINS 

Program BINS is used to smooth the data and to decrease the size of 
th_* data file used for the inversion. The full 512 points analyzed by the 
program TANSMM are far more than needed to describe the occultation adequat- 
ely. Experience has shown that 150 po’ s are sufficient for that purpose. 
Also the inversion program for convenie.. a, always assumes that the occulta- 
tion profile is presented in numerical order in the tangent height. For the 
raw data this ordering hold for the dawn cases but the order must be reversed 
for the dusk cases. 

The user is prompted for the file name in the usual manner. The 
".TH" File is opened and the latitude and longitude are read and printed. The 
user is then prompted for the type of occultation, dawn or dusk. The tangent 
height profile is then > i-ad in the appropriate order from the ".TH" file and 
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the file is closed. The ".OC" file is then opened to read the actual normal- 
ized data counts that correspond to the tangent height profile just read from 
the ".TH" file. 

The program then prompts for the number of points per bin desired. 
For most cases we have found that three is the proper number to choose since 
this allows an in-memeory matrix inversion to be performed for the profile 
inversion.. Since the occultation profile may occur anywhere within the 512 
points, it is desirable to shorten the profile to simplify the data handling. 
The program prompts the user for the number of points to be skipped at the 
beginning of the 512 point block of data and for the number at the end of the 
block. The user may use the printout of the TANSW program as an aid in decid- 
ing vJiich points contain the actual profile. Care must be taken because the 
dusk profiles have been reversed as described above. 

After closing the ".OC" file, a new output file is created with the 
extension ".BN". This contains the binned data points and the corresponding 
tangent heights averaged over the number of binning points (usually three). 
The file is. closed before the program ends. 

6. PROGRAM REALOZ 

This program performs the actual inversion of the occultation profile 
to obtain the profile of the ozone as observed by the UVSP. Using the file- 
name requested, the program opens the corresponding ".BN" file and reads its 
contents before closing it. It then prompts for the ozone absorption cross 
section for the wavelength of the experiment. The geometry matrix is initial- 
ized to zero and other geometric quantities are defined. The main loop to cal- 
culate the geometry matrix is then entered and the triangular matrix is comput- 
ed in the linear array "DELS1 (2601)" which corresponds to a matrix of any size 
up to 51 by 51. The matrix is computed as an equivalent linear array to allow 
the dimensions to be changed easily. The matrix inversion routine MINV is 
then called, the inverse is computed and returned in the same array DELS1. 
The required matrix multiplication is performed by calling the routine GMPRD 
and the ozone profile is printed before the program ends. 
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7. DOUBLE PRECISION FUNCTION DOT 

This function computes the dot product of two Cartesian vectors and 
returns the value of the product. It is called by subroutine ROLLIT. 

t 

8. SUBROUTINE ELLIPS 

This subroutine is called by the program TANSMM to compute the 
height, latitude and longitude of the calculated tangent point which is de- 
scribed as a Cartesian three- vector. The elliptic figure of the earth is used 
to assure the correct altitude. 

9. SUBROUTINE EPHTST 

This routine provides the spline fit of the satellite ephemeris. It 
is called by the program TANSMM. The first step is to initialize the various 
arrays to zero. The ephemeris data is then read from the file (.PT) which was 
opened in TANSMM. The subroutine SPLINE is then called for each of the three 
Cartesian coordinates of the satellite position vector. This provides a cubic 
spline fit which will -be used in function SEVAL as called in TANSMM. Finally, 
the file is closed. 

10. Subroutine GMPRD 

This routine is adapted from the IBM Scientific Subroutine Package. 
It provides the generalized capability to multiply any two compatible mat- 
rices. It is called from the program REALOZ to multiply the inverted geometry 
matrix by the derived observation column vector to obtain the ozone profile. 
As the routine is heavily commented, its operation should be quite clear. 

11. SUBROUTINE LATLON 

This routine calculates the radius, latitude and longitude of any 
point given its Cartesian coordinates. 

12. SUBROUTINE MINV 

This IBM Scientific Subroutine Package routine provides the inverse 
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of any square matrix using the Gauss- Jordan method. The matrix is to be 
passed in the calling sequence as an equivalent linear array to enable chang- 
ing dimensions easily. The comments in the routine should be sufficient for 
rapidly following the structure of the routine. It is cal lied by the program 
REALOZ. 

13. FUNCTION RMXKN 

The purpose of RMXMN is to compute either the maximum of an array of 
values or the minimum of the same array. The choice is made by the value of 
the third parameter in the calling sequence. The method is straight forward 
and uses a brute- force sequence. It should not be used for a very large sized 
array. 

14. SUBROUTINES ROLL IT 

This routine is called from the program TANSMM and provides the rota- 
tion matrices required for the transformation of the various coordinate sys- 
tems. Three matrices are computed. The first, "Rl", moves the z-axis to the 
radius through the satellite. The calling program then performs the transla- 
tion such that the origin is then at the satellite itself. The second matrix, 
"R2", provides the z-axis along the sight- path of the center of the UVSP slit. 
The third matrix, "R3", is used to position the plane of viewing with respect 
to the solar axis of rotation which is used by the sensor as a convenient axis 
of reference. 

15. SUBROUTINE ROTMTX 

This subroutine performs the transformation of a given vector by a 
given rotation matrix. It takes advantage of the properties of orthogonality 
of rotation matrices to enable the inverse transformation using the transpose 
of the matrix as passed in the calling sequence. Which type of tranformation 
is desired, is determined by the value of a flag passed through the calling 
sequence. 
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16. FUNCTION SEVAL 

This double precision function evaluates the cubic spline interpola- 
tion using the spline coefficients computed with subroutine SPLINE. 

17. SUBROUTINE SPLINE 

This subroutine computes the cubic spline coefficients required for 

the interpolation of the satellite ephemeris. Its operation should be clear 
from the comments. 

18. SUBROUTINE SOLE PH 

The solar ephemeris for the year 1980 is computed by evaluating 
several Chebyshev polynomial fits published by the U.S. Naval Observatory. 
Note that it would be necessary to change the coefficients if one were to 

analyze data from another year. This would be a simple matter to do. One 

only has to obtain the coefficients from the Naval Observatory publication 
corresponding to the year of the data. 

19. SUBROUTINE TANHGT 

The Cartesian coordinates of the tangent point for a given line-of- 
sight are computed in double precision by this subroutine. The unit vector of 
the line-of-sight is passed in the calling sequence as well as the distance to 
the tangent point ("TPARAM"). The tangent point coordinates ("XTAN") are re- 
turned together with the single precision magnitude of the radius to the 

tangent point ("RTAN"). 

20. SUBROUTINE UNIVEC 

This routine computes the unit vector from the satellite to the sun 
and returns it in the array ("EUNIT"). 
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c 

c 

c 

c 

c 

c 

c 


1 

c 


c 

c 


1 


c 


100 

c 

15 

exit’ 

c 


c 


c 


c 

99 


program datfnd 

this program reads the “-final data" -files *.fd 
in order to -find by hand the blocks in which the 
occultation occurs, the number o-f the block should 
be noted tor use in program datpgm 

the -following 2 statements must be changed -for Fortran IV 

chararter *30 -fname 

character*6 name 

integc^*2 ismm(256) 

conti nue 

user selects the file name to be examined 
print*, ' enter file specif ication’ 
read*, name 
print*, name 

filename is concatenated with ".fd" extension 
this is unix way of doing it change for rt-11 
f name=name/ / ’ . f d ’ 
pri nt*, fname 

open (9, status=’ ol d ’, f or m= ’ unformatted’ , access=’ direct ’ , 
reel =512, fi 1 e=f name) 

user enters the block number he wants to look at 
print*, ’block number’ 
read*, iblock 

the block is read and written to the screen of the ert 
read (9, rec=iblock) iemm 
print 100, ismm 
format ( lOi 8) 

user selects another block or case or stops 

print*, ’ type 1 for another block, 2 another file, 0 

read the selection 
read*, ians 

i f ( i ans. eq. 1 ) go to 2 
i + ( i ans. eq. 0) go to 99 

user selected another file to examine 
i f ( i ans. ne. 2) go to 15 
close (9) 

user selected another case 
go to 1 

user selects to finish 

cl ose (9) 

stop 

end 
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subroutine spl i ne (n, x , y, b, c , d) ' 

integer n 

real*B x (n) , y (n) , b (n) , c (n) , d (n ) 

ccc 

ccc the coeff icients b(i),c(i), and d(i), i=l,2..,n are computed 

ccc for a cubic interpolating spline 

ccc 

ccc s(x)=y(i) + b (i )*(x-x (i ) ) + c (i ) * <x-x (i ) ) %%2 + 

d (i ) * (x-x (i) )**3 

ccc 

ccc for x(i> .le. x .le. x(i + l) 
ccc 

ccc input 
ccc 

ccc n = the number of data points or knots ( (n. ge. 2) 

ccc x = the abscissas of the knots in strictly increasing order 

ccc y = the ordinates of the knots 

ccc 

ccc output. . 
ccc 

ccc b.c.d = arrays of spline coefficients as defined above, 
cc 

ccc using p to denote differentiation 

ccc - 

ccc y (i ) = s (x (i ) ) 

ccc b (i ) = sp (x (i ) ) 

ccc c(i) = spp(x(i))/2 

ccc d(i> = sppp(x(i))/6 (derivative from the right) 
ccc 

ccc the accompanying function subprogram seval can be used 

cc to evaluate the spline 

ccc 

ccc 

integer nml,ib,i 
real*B t 

ccc 

nml = n-1 

if (n .It. 2) return 
if (n .It. 3) go to 50 

ccc 

ccc set up tridiagonal system 
ccc 

ccc b = diagonal, d = offdiagonal, c = right hand side 
cc 

d(l) = x (2) - x (1) 

c (2) = (y (2) - y ( 1 ) ) /d ( 1 ) 

do 10 i=2,nml 

d ( l ) = x(i + l)-x(i) 

b (i ) = 2.0d00*(d(i-i) + d(i>) 

c(i + l) = (y (i + 1 ) -y (i ) ) /d (i ) 

c(i) = c(i+l)-c(i) 

10 continue 
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t ccc 

ccc end conditions. third derivatives at :a C i > and x (n) 

cc obtained -from divided dif f erences 

ccc 

b C 1 ) = -d ( 1 ) 

b(n) = -d(n-l) 

< c(l) = O.OdOO 

c (n) = O.OdOO 

if (n .eq.3) go to 15 

c(l) = c(3)/(x(4)-x(2))-c(2)/(x<3)— x(l)) 
cm) - c(n-l)/(x (n)-x (n-2) ) - c (n-2) / (x (n-1 ) -x <n-3> ) 
c(l) = c<l)*d(l)**2/(x <4)-x (1)> 
f c (n) = -c(n)*d(n-l)**2/(x (n)-x (n-3)) 

ccc 

ccc forward elimimati'on 
ccc 

15 do 20 i =2, n 

t = d (i-1) /b ( i -1 ) 

$ b (i ) = b (i ) - t*d (i-1 ) 

c (i ) = c (i ) - t*c (i-1 ) 

. 20 continue 

ccc 

ccc back substitution 
ccc 

' c (n) = c In) /b (n) 

do 30 ib=l,nml 
i = n-ib 

c ( ■ ) = (c ( i ) - d(i)*c(i + l))/b(i) 

30 continue 

( CCC 

ccc c(i) is now the sigma (i) of the text 
ccc 

ccc compute polynomial coefficients 
ccc 

b(n) = (y(n) - y (nml ) ) /d (nml ) + d (nml ) * (c (nml ) + 
2.0d00*c(n)) 

do 40 i=l,nml 

b ( i ) = (y (i +1) - y(i))/d(i) - d<i)*(c(i+l) + 2.*c(i)) 
d ( i ) = (c (i +1 ) - c(i))/d<i) 
c (i ) = 3.0d00*c(i) 

40 continue 

( c (n) = 3.0d00*c(n) 

d (n) = d (n-1 ) 
return 

ccc 

50 b (1) = < y ( 2) — y ( 1 ) ) / ( x ( 2 ) -x ( 1 ) ) . 
c ( 1 ) = O.OdOO 
d ( 1) = O.OdOO 
b (2) = b ( 1 ) 
c (2) •-= O.OdOO 
d (2) =0. OdOO 
return 
end 

< 22 ) 
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■function rmxmn (nd, data, ixn) 
dimension data(l) 

c returns max value (ixn=l) or min value <ixn=0) of array 

rmxmn=data ( 1 ) 
if (ixn.eq.O) go to 200 
do 100 j=2,nd 

if (data (j) . gt. rmxmn) rmxmn=data ( j ) 

100 continue 

return 

200 do 210 j=2, nd 

i f (data ( j ) . 1 1. rmxmn) rmxmn=data ( j) 

210 continue 

return 
end 
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subroutine uni vec (xsat, xsun, eunit) 
returns double precision unit vector 
along the vector (xsun-xsat) 
real *8 xsat (3) , xsun (3) , eunit (3) , xmag 
xmag=0. OdOO 
do 10 i-1,3 

xmag=xmag-*- (xsun (i ) -xsat (i > ) **2 
continue 
xmag=dsqrt (xmag) 
do 20 1*1,3 

euni t (i ) = (xsun (i ) -xsat (i ) ) /xmag 

continue 

return 
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subroutine tanhgt (:<sat , euni t , tparam, xtan, rtan) 
c returns the location oi the tangent point 

real^B d, xsat (3) , euni t (3) , tparam, xtan (3) 
real *8 rtan 

d= ( (xsat (2) fteuni t (3) -xsat (3) fteuni t (2) ) **2+ 

1 (xsat (3) *euni t ( 1) -x sat ( 1 ) fceuni t (3) ) »*2+ 

2 (xsat ( 1 ) ieuni t (2) -xsat (2) leuni t ( 1 ) ) * * 2 ) 
tparam=dsqrt (xsat ( 1 ) #*2+xsat (2) #*2+::sat (3) **2-d) 
b=dsqrt (d) 

do 10 i=l,3 

xtan (i ) =xsat (i ) +tparam*euni t (i ) 

10 continue 

rtan=sngl (d) 

return 

end 
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subroutine cross (x , y, xy) 

real *8 :< (3) , y (3) , xy (3) 

xy(l) = x(2)*y<3) - x(3)*y 2) 

xy (2) = x (3) *y ( 1 ) - x(l)ty(3> 

x y ( 3 ) = x < 1 ) *y (2) - x(2)*y(l) 

return 

end 
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program tansmm 

r 

c program to compute the tangent points For 

c the occultation blocks 

c 

common/xyz/tt (50) , xx (50) , yy (50) ,zz (50) , xcl (50) , ::c2 (50) , xc3 (50) , 
xycl (50) ,yc2<50) ,yc3(50) ,zcl (50) , zc2 (50) , zc3 (50) 
real *Btt, xx , yy, zz, xcl, xc2, xc3, ycl , yc2, yc3, zc 1 , zc2, zc3 
i nteger *2i month, i out, i check, i dayyr 
c the Following statement is a UNIX space saver 

c it may be replaced by a simple dimension statement 

c For use on other systems 

c 

automatic htan(512) 

real *Bt star t , tbl ock , t , xsun (3) , xsat (3) , euni t (3) ,xtan (3) , tparam, 
xti nc, srad, si at , si on, pitch, yaw, roll, xdum (3) , ydum (3) , r 1 (3, 3) , 

:: r2<3,3) 

real *0sol axs (3) , aphi , bphi , tphi , eel ptc, r3 (3, 3) , sol tl t , radg, sundec , 
xpi , sunzen , rtan , tanl at , tanl on 
character*30 Fname 
character *6 name 
character*6 Fxtim(7) 

c the Following experiment numLers need special oFFsets 

c thus they are singled out in a data statement 

data Fxtim/’ VI 1134', ’ VI 11 40’ , ’VI 1144’ , ’ VI 1 171 ’ , ' VI 1 172’ , 

1 '011173’ , 'VI 1174’ / 

data eel ptc/23. 4333333333d00/ 
data sol tl t/7. 25d00/ 
data aphi /73. 666666666d00/ 
data bphi /l. 3950333333d -02/ 
radg = datan ( 1 . dOO) /45.d00 
pi = lBO.dOOtradg 
23000 continue 

c user supplies the File or experiment name 

print*, ' enter File speciFication' 
read *, name 
print*, name 

c Following statement concatenates extension to Filename 

Fname=name//' .oc' 
print*, Fname 

c open the File with the occultation data 

open (9, status*' ol d ' , Fi 1 e*Fname, access*' sequential ' , 

1 reel =00, Form*’ Formatted' ) 

c File needs to be placed at beginning oF information 

rewind 9 

c read the date and other time values 

r^ad (9, 10) idayyr, imonth 

read (9, 20) tstart , tbl ock, pi tch, yaw, rol 1 , tine 
close (9) 

c concatenate to get the partial ephemeris data 

Fname = name//".pt" 
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open (9, -form ■ "formatted*, status » "old", -file ■ fname) 
rewind 9 

c get the spline fit to the ephemeris 

call ephtst (tblock, npt) 
soltlt = soltltftradg 
eclptc * eclptcfcradg 

tphi = dbl e (1 980-1B50) + (dbl e ( i dayyr ) +tbi ock/24d00) /365. 24d00 
tphi ■ (aphi +bphi *tphi ) *radg 
c compute the solar rotation axis 

solaxs(l) = dsin (soltlt) *dcos (tphi ) 

solaxs(2) = dsin (sol tl t) *dsin (tphi ) *dcos (eclptc) -dcos (sol tl t) * 
xdsin (eclptc) 

solax«(3) « dsi n (sol tl t ) *dsi n (tphi ) tdsin (eel ptc) +dcos (sol tl t) * 
xdcos (eclptc) 
iout * 1 
i check = 0 
t - tblock 

c compute the special offset if necessary 

do 200 ii«*l,7 

if (name .eq. fxtim(ii)) t“t-262. 144/3. 6c03 
200 continue 

i-dayl = i dayyr 

c loop over, the two occultation blocks (512 points) 

do 23003 i“= 1,512 
i dayyr=i day 1 

c get the solar ephemeris for the solar location 

call sol eph ( t , i dayyr, :<sun, sun cec ) 
sundec = sundeeftradg 

c use the spline fit to get the satellite location 

c by calls to the function seval 

xsat(l) = seval (npt, t, tt, xx , xc 1 , xc2, xc3) 
xsat(2) *= seval (npt , t, tt, yy, yc 1 , yc2, yc3) 
xsat (3) = seval (npt , t., tt, zz , zc 1 , zc2, zc3) 
c call latlon to get radius, latitude & longitude 

cal 1 1 at 1 on (xsat, t , srad, si at, si on) 
c get the unit vector to the sun center 

call uni vec (xsat, xsun, euni t) 
c call the rotation matrix forming routine 

call rollit (xsat , srad, si at , si on, pi tch, yaw, roll, euni t, 
x sol ax s, r 1 , r2, r3, sunzen) 

101 format (3(2x, lpgl4. 7) ) 

c call the tangent height calculator routine 

call tanhgt (xsat, euni t , t par am, xtan, rtan) 
c perform the rotation transf ormati on 

call rotmtx (xsat, r 1 , xdum, 1 ) 
c do i. translation 

xdum (3) = xdum (3) -srad 
c oo another rotation 

call rotmtx (xdum, r 2, ydum, 1 ) 
c compute the tangent point coordinates 

xtan(3) = tparam 

xtan(2) = -tparam*dsi n (yaw) /dcos (yaw) 

28 ) 
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xtan(l) = -tparam*dsi n (pi tch) /dcos (pi tch) 
c perform two more rotations 

call rotmtx (xtan, r3, ydum, -1 ) 
call rotmt:: (ydum, r2, xdum, -1 ) 
xdum ( 1 ) = -xdum(l) 
xdum(2) = -xdum(2) 
xdum (3) = srad+xdum (3) 
c do another rotation 

call rotmtx (xdum, rl, xtan, -1) 
c compute the unit vector 

call uni vec ( x sat , t an , eun i t ) 
c get the tangent point 

call tanhgt (:<sat, euni t , tparam, xtan, r tan) 
c compute position on elliptical earth 

cal 1 el 1 ips (xtan, t, rtan, h tan (i ) , tanl at , tan Ion, sundec, sunzen, 
xpi ) 

c print, every tenth point 

if (.not. (mod ( (i+9) , 10) .eq.O) )goto 23005 
pri nt * , "htan (i ) = “ , i , htan ( i ) , tanl at , tanl on 
23005 continue 

t = t+tinc 
iout = iout+1 
23003 continue 

c open the "tangent height -file * . th 

fname = name//".th" 

open (9, status = "new" -ile = -fname, access = "sequential ", reel 
xBl,form = "-formatted") 

c write the tangent point coordinate on *.th 

wr i t e ( 9 , 1 0 ) i ou t 
wr: te (9, 30) tanl at , tanl on 
wri te (9, 40) htan 
cl ose (9) 
stop 


10 

■format (1615) 


20 

format (4d20. 

12) 

30 

f ormat (Bf 10 

.3) 

40 

format (8f 10 

.5) 
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program bins 
c 

c this program performs a binning and averaging 

c of the occultation being considered 

c 

integer *2 

i out, j, j j , nptbin, nski p, i start, iend, i , nendsk, i pick- 
real *B tstart, tblock, pi tch, yaw, rol 1 , tine 
real *4 htan (512) , data (512) , hbi n (51 ) , dbi n (51) 
character*3C» fname 
character *6 name 

ccc specify and open input file *.th 

pri nt*, ’ enter file specification’ 
read*, name 
f name=name//’ . th’ 


open (9, 
1 


ccc 


c 


01 

c 

c 


status=’ old’ , ecce55= J sequential ’ , form=’ formatted’ , 
reel =S1 , f i 1 e=f name) 
rewind 9 

read input and reverse array if necessary 
read (9, 10) iout 
read(9»30> tanlat, tanl on 
print*,’! out =’ , i out , tanlat, tanl on 
user inputs whether dawn or dusk case 
print*,’ dawn or dusk, oawn=l , dusk=2’ 
read*, i pi ck 

if (ipick. eq. 1 ) read (9, 40, end=201 ) (htan (i ) , i =1 , i out) 
if (ipick. eq. 2) reaa (9, 40, er.d=201 ) (htan (i ) , i =i out ,1,-1) 
ci ose (9) 

read the file containing the occultation data named 
* . oc 


f name=name//’ . oc’ 
print*, fname 

open (B, status=’ old’, access=’ sequential ’ , f orm=’ f or matted 
1 reel =80, f i le=f name) 
rewind B 

read (8, 10) i dayyr, imonth 

read (B, 20) tstart, tbl ock, pi tch, yaw, roll, tine 
i f (ipi ck. eq . 1 ) read (8, 41 , end=202) (data ( j ), j=l , i out ) 
i f ( ipi ck. eq. 2) read (8, 41 , end=202) (data ( j ), j=i out, 1 , -1 ) 
02 close (B) 

ccc s c up bins 

c user supplies the number of points per bin 

50 print*, ’enter number of paints per bin’ 
read*, nptbin 

c user supplies the number of points to skip 

c at the beginning of the raw data file *.oc 

print*, ’ enter number of points to skip at beginning’ 
read*, nski p 

user supplies the number of points to skip at 
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c the end of the raw data -file 

pri nt*, ’ enter number of points to skip at end’ 
read*, nendsk 
i out=i out-nendsk 
i start=nski p+1 

c prepare the file which is to contain the binned 

c occultation data in a file *.bn 

fname=name//’ .bn’ 

□pen (7, status^’ new’ , access=' sequential * , for m=’ formatted ' , 
1 reel =80, file=f name) 

write (7, 42) tanl at , tanl on 
i end=i start+npt bi n-1 
j j = l 

55 hbin(jj)=0. 0 
dbin ( j j) =0. 0 

c compute the averages of the tangent heights and data 

do 60 i =i start, i end 


hbi r. ( j j ) =hbi n ( j j ) +htan (i ) 
dbi n ( j j ) =dbi r. ( j j ) +bata (i ) 

60 continue 

hbi n ( j j ) =hbi n ( j j ) /r.ptbi n 

dbi n ( j j ) =dbi n ( j j ) /nptbi n 

write(7,43) j j , hbi n ( j j ) , dbi n ( j j ) 

i start=i start+nptbin 

iend=iend+nptbin 

i f ( i end. gt. i out ) go to 99 

jj=jj + l 

go to 55 

99 continue 

ciose(7) 
stop 

10 format(16i5) 

20 format (4d20. 12) 

30 format (Bf 10 .3) 

40 format (BflO .5) 

41 f ormat (Bf 10. 7) 

42 f ormat (2f 10. 3) 

43 f ormat (i5,f 10.3, f 12.7) 
end 
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c 

c 

c 

c 

c 

c 

c 


ccc 

1 


c 

c 

c 


c 


1 

c 


c 

1 

c 


c 

1 


c 

c 


1 


c 

c 


program aotopt 

this program reads the SMM ephemeris data -file 
and pulls out the part of the ephemeris which 
corresponds to the occultation period 
this is placed on a file with extension 7 .pt’ 
the input comes from a file with extension ’.ao’ 

real *8 exptno, two, z date (2) , zsec (2) , zday (2) , zz (56) 
real *8 xdate, xday, xsec, xdut 

real *B x ( 10) , y ( 10) , z ( 10) , vx ( 10) , vy ( 10) , vz (10) 

character *30 fname,pname 

character *6 name 

specify and open input file 

continue 

pri nt*, 7 enter file specification 7 

read *,name 

if (name .eq. 7 z 7 ) stop 

the next statement concatenates name with the . ao 
extension to qet the ephemeris file required, 
f name=name//’ .ao’ 

concatenate for the shortened data file *.pt 
pname=name// 7 . pt 7 
print *,fname, pname 
open the epehmeris file for reading 
open (9, status= 7 ol d 7 , access= 7 direct 7 , 
reel =512, f i 1 e=f name) 

get the times of interest from the input. 
print*, 7 input start and end time in sec 7 
react, tstart , tend 

read *.ao to get experiment number etc. 
read (9, rec=l ) exptno, two, 

(zdate (i ) , zday ( i ) , zsec (i ) , i =1 , 2) , zz 
write it to standard output 
print 101, exptno, (zdate (i ), zday (i ), zsec (i ), i =1 , 2) 
i =1 
i j k = 1 

open the shortened ephemeris file for output 
open (B, ststus= 7 new 7 , reel =252, f i 1 e=pname, 
access= 7 sequential 7 , f orm= 7 formatted 7 ) 
start the read loop 
continue 
i =i +1 

read the ephemeris data 

and print the data it time between tstart and tend 
read (9, rec=i , end=90) xdate, xday, xsec, xdut 

, (x (k) ,y(k) ,z (k) , vx (k),vy(k),vz (k) ,k=l, 10) 
i f ( xsec. 1 t. tstart . or . xsec. gt. tend) go to 2 
printl02, xdate, xday, xsec, xdut 
write out every tenth data point 
as epehemeris is needed every ten 
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c since the spline -fit will be sufficient for 

c finer time detail 

k=l 

print 105, xsec, x (k) , y (k) , z ( k ) 
write (8,105) xsec, x (k),y(k),z(k) 

105 format (f 6. 0, 3 (2x , dl7. 10) ) 

xsec=xsec+2. 0 
3 continue 

print 104, ijk 
104 format(ilO) 

i jk=i jk+1 

if(ijk .ge. 51) goto 10 
goto2 

90 print*, ' eof ' 

10 continue 

cl ose (9) 
cl ose (8) 
gotol 

101 format (20x, 'experiment number ',f9.1/, 

1 ' start date = ' , f 9. 1 , 2x , f 5. 1 , ' th day of year at 
’ , f 9. 1 , 

2 .’secs of day'/' end date = ' , f 9. 1 , 2x , f 5. 1 , ’ th day of 
year ' 

3 'at ' , f 9^ 1 / ) 

102 format(/' date= ',f9. 1,' day= ',f5.1,' secs= ',f9.1, 

1 ’ dut= ’ ,f 14.5/10x, lhx. 19x, lhy, 19x, lhz/9x,2hvx, lBx,2hvy, 

2 lBx , 2hvz ) 

103 f ormat (i3,3f20. 13/3x3f 20. 13) 
end 
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subroutine minv( a,n,d,l,m ) 
c 


c 

c subroutine minv 

c 

c purpose 

c invert a matrix 

c 

c usage 

c call minv (a, n, d, 1 , m> 

c 

c description of parameters 

c a - input matrix, destroyed in computation and 

replaced by 

c resultant inverse, 

c n - order of matrix a 

c d - resultant determinant 

c 1 - work vector of length n 

c m - work vector of length n 

c 

c remarks = 

c matrix a must be a general matrix 

c 

c subroutines and function subprograms required 

none 


c 

c method 

c the standard gauss-jordan method is used, the 

determinant 

c is also calculated, a determinant of zero indicates 

that 

c the matrix is singular. 

c 

c 


c 

c subroutine mi nv (a, n, d, 1 , m) 

dimension a ( 1 ) , 1 ( 1 ) , m ( 1 ) 
c 
c 


c 

c if a double precision version of this routine is 

desired, the 

c c in column 1 should be removed from the double 

preci sion 

c statement which follows, 

c 

c double precision a, d, bi ga, hoi d, dabs 
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c 

c the c must also be removed -from double precision 

statements 

c appearing in other routines used in conjunction with 

thi s 

c routine, 

c 

c the double precision version o-f this subroutine must 

al so 

c contain double precision -fortran -functions. abs in 

statement 

c 10 must be changed to dabs. 

c 

c 


c 

c search -for largest element 

c 

o=i.o 

nk=-n 

do BO k = i,n 
nk.=nk+n 
1 ( k ) =k 
m Ck) =k ~ 
kk=r k+k 
bi ga=a (kk) 
do 20 j=k,n 
i z=n* ( j — 1 ) 
do 20 i=k,n 
i j=iz+i 

10 i -f ( abs(biga)- abs(a(ij>)) 15,20,20 
15 biga=a(ij) 

I (k)=i 
m ( k ) = j 
20 continue 
c 

c interchange rows 

c 


j=2 (k) 

i -f ( j-k) 35,35,25 
ki =k-n 
do 30 i=l,n 
ki =ki +n 
hold=-a (ki ) 


ji =ki -k+j 
a ( ki ) =a ( ji ) 

30 a( ji) =hol d 
c 

c interchange columns 

c 

35 i=m(k) 

if(i-k) 45, 45, 3B 


\ 
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3B jp=n*(i-l) 
do 40 j=l,n 
jk=nk+ j 
ji=JP + j 
hoi d=-a ( jk) 
a ( jk) =a ( ji ) 

40 a ( ji ) =hold 

c 

c divide column by minus pivot (value of pivot element is 

c contained in biga) 

c 

45 if (biga) 48,46,48 
46 d=0. 0 

print*, ’d=0’,k, biga 
return 

48 do 55 i=l,n 

if(i-k) 50,55,50 
50 ik=nk+i 

a<ik)=a(ik)/(-biga) 

55 continue 
c 

c • reduce matrix 

c 

do 65 i=l,n 
i k=nk+i 
hoi d=a ( i k) 
i j=i -n 
do 65 j=i,n 
i j=i j+n 

if(i-k) 60,65,60 
60 i-f (j-k) 62,65,62 
62 kj=ij-i+k 

a(ij)=hold*a(kj)+a(ij) 

65 continue 
c 

c divide row by pivot 

c 

k j=k-n 
do 75 j=l,n 
k j=k j+n 

if ( j-k) 70,75,70 
70 a (k j ) =a (k j ) /bi ga 
75 continue 
c 

c product of pivots 

c print *, ’ mi nv’ , k, d , bi ga 

c 

d=d*bi ga 

replace pivot by reciprocal 


c 

c 

c 


a (kk) =1 . 0/bi ga 
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BO continue 

■final row and column interchange 
print*, ’got to r/c inter’ 

k=n 

100 k= (k-1) 

i f (k) 150,150,105 
105 i =1 (k) 

if(i-k) 120,120,108 
108 jq=n* (k-1 ) 
jr=n* (i -1 ) 
do 110 j=l,n 
jk=jq+j 
hoi d=a ( jk) 
j i = jr+ j 
a ( jk) =-a ( ji ) 

110 a ( ji ) =hol d 
120 j=m(k) 

if ( j-k) 100,100,125 
125 ki =k-n 

do 130 i=l,n 
ki=ki+n _ 
hoi d=a (ki ) 
ji =ki -k+ j 
a (ki ) =-a ( ji ) 

130 a ( ji ) =hol d 
go to 100 
150 return 
end 
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subroutine ephtst (tblock,npt) 
c 

c this routir ■* reads the ephemeris data and 

c performs a spline -fit for the cartesian coord- 

c ates which will later be evaluated for given times 

c by function seval 

c 

common/xyz/t (50) , xx (50) , yy (50) , zz (50) , xcl (50) , xc2 (50) , 

1 xc3 (50) , ycl (50) , yc2 (50) , yc3 (50) , zcl (50) , zc2 (50) , 

2 zc3 (50) 

real *8 t, xx , yy, zz , xcl , xc2, xc3, ycl , yc2, yc3, zcl , zc2, zc3 
real *4 tt(50) 

realiQ tblock 

c do loop to zero the arrays as initial values 

do 5 i = 1 , 50 
tt (i )=00.00 
xx ( i ) =00. dOO 
y y ( i ) =x x ( i ) 
z z ( i ) =x x ( i ) 
t (i ) =xx ( i ) 

•xcl ( i ) =x x (i) 
xc2 (i ) =xx_(i ) 
xc3 (i ) =xx (i ) 
ycl (i ) =x:< (i ) 
yc2 (i ) =xx (i ) 
yc3 (i ) =xx (i ) 
zcl (i)=xx (i> 
zc2 (i ) =xx (i ) 
zc3 (i ) =xx (i ) 

5 continue 

npt=l 

c read the ephemeris data to prepare for the spline fit 

do 100 k=l , 50 

read (9, 2000, end=300) tt(k),xx(k),yy(k),zz(k) 

2000 format (f 6. 0, 2x , dl7. 10, 2x, dl7. 10, 2x , dl7. 10) 
npt=npt+l 

c convert units to current problem 

t (k) =dbl e (tt (k) ) /3600. dOO 
xx <k)=xx (k)*l.d4 
yy ( k) =yy (k) *1 . d4 
zz (k)=zz ( k) * 1 . d 4 
100 continue 
300 npt=npt-l 

c call the spline routine three' times 

c to fit each of the cartesian coordinates 

call spl i ne (npt , t , xx , xcl , xc2, xc3) 
cal 1 spline(npt,t,yy,ycl,yc2,yc3) 
cal 1 spline (npt, t,zz,zcl,zc2,zc3) 
c print the values used 

print*, ’npt=’, npt 
c 
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da BOO kkk=l,npt 
SOI -f or mat (4g IB. 1 1 ) 

print BOl , t (kkk) , s<j< (kkk) , yy (kkk) , zz (kkk) 
802 -format (lBx,3glB. 11) 

BOO continue 
close (9) 
return 
end 
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pro?'’' am realoz 

this program inverts the geometry matrix 
to obtain the ozone pro-file 

dimensi on ans (51 ) , glamda (51) , tau (51 ) , 1 work (51) , mwork (51 ) 

1 , htan (51 ) , havg (51 ) 

dimension oelsl(2601) 
character#30 f name 
character*6 name 

arithmetic statement -function for the figure 
of the earth 

geoid (x ) *=6378.388# ( 1 . -. 3367003e-2#si n (x ) **2+. 70B5e-5* 
lsin (:<*2) *#2> 

rdd=3. 14159/1B0. 
continue 

pri nt*, ' enter file specification' 
read#, name 

concatenate for the binned data file 
f name=name//’ .bn’ 

• open ( 9, status*’ old', access*’ sequent; al ' , f orm*’ formatted' , 
1 recl*=B0, f i 1 e=f name) 

read (9, 42) tanl at , tanl on 
format (2f 10.3) 
do 57 j k q = 1 , 5 1 

read (9, 43, end=5B) jjj,htan(jjj),tau(jjj) 

format (i 5, f 10.3, f 12.7) 

conti nue 

close (9) 

ntan=jjj 

print*, tau 

user supplies ozone cross section 
print*,' input cross section... ’ 
read *,csx 

print #, ’ cross section used = ',csx 

zero out the matrix initially 
do 5 j=l , 2601 
del si ( j ) *=0. 

hupper s, htan (ntan) 
dhtan=htan (ntan) -htan (ntan-1 ) 
ntar=ntan -1 
re=geoi d (tanl at*rdd) 
r upper =re+hupper+dh tan 
ntan2*=ntan#ntan 
i =0 

j*=0 

continue 0 

j=i *ntan+i +1 
i *=i -*-1 

rtan*htan (i ) +re 
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havg (i ) ** (htan (i ) +htan (i +1 ) ) /2. 

h“htan (i ) 

rtan2*rtan*r tan 

rupper^re+huppar-'-dhtan 

start=0. 0 

k®i -1 

j j“j 

9 k=k+l 

r=re+htan (k+1 ) 

sdi st=sqrt (r#r — rtan2) 

del si <jj)“ (sdi st-start ) *1 . 0e-02 

start=sdi st 

j j«ntan+ i j 

It. ntan2) go to 9 
glamda (i ) =0. 0 

i( (tau(i) .It. l.e-36) go to 15 
gl amca (i ) =0. Stalog ( 1 . /tan (i ) ) /csx 

15 continue 

if(j .It. ntan2+l) go to B 

20 continue 

21 continue 

c . call the matrix inversion routine 

call minv (del si , ntan, d, 1 work, mwork) 

c ‘ call thd matrix multiplication routine 

call gmprd (del si . gl amoa, ans, ntan , ntan, 1 ) 
print 991 , tanl at, tar.l on 

991 4 ormat (2x , ” tangent area latitude= ’,46.2, ’ long.*= 

’,47. 2) 

do2139 nn=l,ntan 

2139 ans (nn) =abs (ans(nn))*l. 0e-07 
do 21 AO nn=*l,ntan,3 

pri nt992, 

havg (nn) , ans (nn) , havg (nn+1 ) , ans (nn+1 ) , havg (r.n+2) . 

1 ans(nn+2) 

992 -format (3 (Opt 6. 2, 2x , IpelO. 2, 4x ) ) 

2140 continue 
stop 
end 
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subroutine gmprd (a, b, r, n, m, 1 ) 


subroutine gmprd 


purpose 

niultiply two general matrice 


s to form a resultant 


matri x 


usage 

call gmprd (a, b, r , n, m, 1 ) 


description of parameters 

a - name of first input matrix 
b - name of second input matrix 
r - name of outout matrix 
n - number of rows in a 

m - number of columns in a and rows in b 
1 - number of columns in b 


remarks ^ 

all matrices must be stored as general matrices 
matrix r cannot be in the same location as matrix a 
matrix r cannot be in the same location as matrix b 
number of columns of matrix a must be equal to number 


of matrix b 


subroutines and function subprograms required 


method 

the m by 1 matrix b is premultiplied by the n by m 


and the result is stored in the n by 1 matrix r. 


subroutine gmprd (a, b , r , n, m, 1 ) 
dimension a(l),b(l),r(l) 


ir-0 
i k«-m 

do 10 k = l,l 
i k=i k+m 
do 10 j c l,n 
ir=ir+l 
ji-j-n 


« I 
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r C i r ) =0 
do 10 i=l,m 
ji=ji+n 
i b=ib+l 

10 r (ir) =r (ir) +a ( ji ) *b (i b) 
return 
end 


c 

c 

c 

c 

c 


1 

2 


1 

2 


c 

c 


ccc 

c 


100 


101 
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program datpgm 


this program reads the occultation 

data and culls out only the part of interest 

which is contained on two DEC blocks 


real*4 strtm, stoptm, data (512) 
real *8 si amst , si amnd , ::1 amst , :<1 amnd, aO, al, a2 
real 48 t start , tbl ock, pi tch, yaw, roll,pi,tinc,tbl sec 
integer *2 i smm (256) , i smml (256) , i smm-f 1 (512) 
equi val ence (i smm ( 1 ) , i smm-f 1 ( 1) ) , (i smml ( 1 ) , i smm-f 1 (257) ) 
integer*2 i 1 , i 2, i 3, i 4, i 5, i 6, i 7, i 13, 122, i 23, i24, 
i 25, i 26, i 27, i 28, i 29, i 80, i 81 , i 100, i 109, i 1 10, i 1 1 1 , 
ill2,ill3,ii 4879 (32) , i i 160 (96) , y 


•filename without 
speci f i cat i on ’ 


ex tens! on 


integer *4 i 2021 , x , i s67, i s69 
integer *4 ismm4(128) 

real *4 i 101 , i 103, i 105, i 107, i 1 14, i 1 16. i 1 1C, i 120, 
i 122, i 124,1 126, i 128, i 130,1 132, i 134, i 136, i 138, i 140, 
i 142, i 144, i 146, i 148, i 150 
character*30 fname 
character 46 name 
continue 

user supplies the 
print*, * enter -file 
read *, name 
pri nt*, name 

open files required for this run 
first the final data file 
fnafne=name//’ .f d’ 
pri nt*, fname 

read block number and file number 
user supplies the header block number 
pri nt*, ’ header block number’ 
read*, i hbl ck 
print*, ’ ibl ock, ipage’ 
read*, i block, ipage 

open (9, status=’old’ , form=’ unformatted ’ , access=’ di rect ! 
reel =512, f i 1 e= fname) 

read the final data file 
read (9, rec=l ) i smm 
read (9, rec=l ) ismm4 
print the data read 

print 100, i smm ( 10) , ismm ( 1 1 ) , i smm ( 12) , i smm (13) 
format(’ start time ’,i3,’ ’ ,-i2, ’ : ’ , i2, ’ : ’ , i3) 
print 101, i smm ( 1 4) , i smm ( 15) , i smm ( 16) , i smm ( 17) 
formatt’ stop time ’,i3,’ ’ , i 2, ’ : ’ , i 2, ’ : ’ , i 3) 
strtm=3600. * ( i smm (ID) +60. * ( i smm ( 12) ) +i smm (13) 
5toptm=3600. * (ismm ( 15) ) +60. * ( i smm (16) ) +i smm (17) 


(2,18, etc) 


print*, 
print*, 
print*, ’ 


start time-secs 
’stop time -sec 
experiment type 


i smm (14) 


i smm (60) 


i smm ( 10) , ’ 

’ ’ , stoptm 


strtm 


D 
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print*,’ no. o-f detectors used ’,'ismm(61) 
print*, i smm4 (34) , i smm4 (35) 
print*, (ismmdsq) ,irq=£5,72) 
c compute the actual wavelengths in angstroms 

xl amst=dbl e (i smm4 (34) ) 
xl amnd=dbl e (i smm4 (35) ) 

print*,’ xlamst =’ , x 1 amst, ’ xlamnd =’, xlamnd 

aO=l . 9707373047d03 
al=5. 128015298d-3 
a2=2. B7B6d— 10 

slamst=2.0* (aO- (xl amst+206. ) *al - ( (x 1 amst+206) **2) *a2) 
si amnd=2. 0* (aO- (xl amnd+206. ) *al- ( (xl amr.d+206. ) **2) *a2) 
pr i nt*, ’ si amst =’ , si amst , ’ si amnd =’,slamnd 
print*,’ no. o-f wld steps ’,ismm(7B) 
print*,’ wld step size (dw) ’.ismm(79) 
read (9, rec=ihbl ck ) 

i 1 , i2, i 3, i 4, i 5, i 6, i 7, x , x , y, i 13, x , x , x , i 2021 , 

1 i 22, i 23, i 24, i 25, i 26, i 27, 1 28, i 29, x , x , x , x , x , x , x , x , x , 

2 i i 4879, i BO, i Bl , x , x , x , x , x , x , x , x , x , i 100, 1 101 , i 103, i 1 05, i 107, 

3 il09, illO, ill, 1112, ill 3, 1114, 1116, i 118,1 120, 

4 i 122, 1 124, 1 126, i 128, 1 130, 1132, 1 134, 11136, 1 13B, 

5 i 1440, i 1 42, 1 1 44, 1 146, i 148, i 150, x , x , x , x , i i 1 60 

1 cayyr=i 29 
1 hr=i 25 = 


i mi n=i . 


i sec=i 27 
1 msec=i 28 

c compute roll pitch yaw etc -from the .-fd file 

pi =4. d00*datar. ( 1 . dOO) 
pitch = dbl e (ismm (23) ) /10. d00/3600. dOO 
yaw = dbl e ( i smm (24) ) / 10. d00/3600. dOO 
roll = dbl e (i smm ( 25) ) /lOO.dOO 
pitch = pi tch*pi /1B0. dOO 
yaw = yaw*pi / 180. dOO 
roll = rol 1 *pi /ISO. dOO 
c compute start and end times 

tstart=dbl e (i hr ) +(dble(imin) *6dl+dbl e ( i sec) ) /3600d0+ 

1 cbl e (i msec ) /3. 6d6 
c get proper offsets 

i f (ihbl ck. eq. 2) i i nc=i bl ock-3 
if (ihblck.eq. 19) i i nc=i bl ock-20 
if (ihblck.eq. 3B) i i nc=i bl ock-3B 
if (ihblck. eq. 53) i i nc=i bl ock-56 
i f ( i hbl ck. ne. 2. and. i hbl ck . r.e. 19) print*, ’ error in header 
block numb 
1 er ’ 

c compute proper block start time 

tbl ock=t start +dbl e ( . 01 6* i 101*256* (i inc) ) /3. 6d3 
t i nc= . 0 1 6dOC»* db 1 e ( i 1 0 1 ) /3600 . dOO 
c print the data for convenience 

print 200, i 100 

200 f ormat ( lx , ’ no. of actual data points found ’,i5) 
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print 201,1101 

201 f ormat ( lx ,’ mean time between data points in inner’ 

1 ’loop in units o-f 16ms’, lpel2.5) 

print 202, i 103 

202 format (lx ,’ max. time gap found ’,lpel2.5) 
print 203, il05 

203 f ormat ( lx , ’ mi n time gap f oui.d’ , l-.,el2.5) 
print 204, il05 

204 format ( lx ,’ standard devi ati on ( ampl e) ’ , lpel2. 5) 
print 205, i 10? 

205 f ormat ( lx , ’ no. of inner loops’,i5) 
print 206, i 110 

206 format ( lx, ’ count of next loop’,i5) 
print 206, ill 1 

print 206, i 112 
print 207,i 113 

207 f ormat (lx , ’ count of executions in this record ’,i 5) 
print 20B , i 1 1 4 

20B format (lx ,’ mean time be uw en inner loops’, lpel2.5) 
print 209, 1116 

209 f ormat ( lx ,’ mean time between 2nd 1 oops’ , lpel2. 5) 
print 210, i 1 IB 

210 f Drmat ( lx 3 ’ mean time between 3rd 1 oops’ , lpel2. 5) 
print 211,il20 

211 f ormat ( lx ,’ mean time between outer 1 oops’ , lpel2. 5) 
print 212, i 122 

212 f ormat ( lx ,’ mean time between executions’ , lpel2. 5) 
print 213, i 124 

213 format (lx, ’ max time between inner 1 oops’ , lpel2. 5) 
print 214,1126 

214 f ormat ( lx, ’ max for next ’ , lpel2. 5) 
print 214, i 128 

print 223, i 130 

223 f ormat ( lx ,’ max for outer ’ , 1 pel2. 5) 

print 215,il32 

215 f ormat ( lx ,’ max for executions’ , lpel2. 5) 
print 216, i 134 

216 f ormat ( lx , ’ mi n time between inner 1 oops’ , lpel2. 5) 
print 217, i 136 

217 f ormat ( lx , ’ mi n for next ’ , lpel2. 5) 
print 217, i 13B 

print 21 B, i 140 

218 f ormat ( lx , ’ mi n for outer’ , lpel2. 5) 
print 219, i 142 

219 f ormat ( lx , ’ mi n for executions’ , lpel2. 5) 
print 220, i 144 

220 format ( lx ,’ standard deviation for inner loop mean 
time’ , lpel2. 5) 

print 221, i 146 

221 format (lx ,’ standard deviation for next loop’ , lpel2. 5) 
print 221 , i 148 

print 222, i 150 

46 
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222 -format ( lx , ’ standard deviation -for outer 1 oop’ , lpel2. 5) 
print*,’ i 25, 6, 7, 8, 9 ’ , i 25, i 26, i 27, i 2B, i 29 

c now read the first block of the occultation 

read (9,rec=iblock) ismm 
c then read the second block 

read (9, rec=i bl ock+1 ) ismml 
close (9) 

c prepare to open the *.oc 

c file for writing the occultation blocks 

fname=nane//’ . oc’ 

c convert to floating point 

do 10 i =1 , 512 
10 data (i ) =i smmf 1 (i > 

c get max value of data array 

datmax = rmxmn (512, data, 1 ) 
c normalise data array to datmax 

do 20 i =1 , 512 

20 oata(i) = data (i ) /datma:: 

open (9, status=’ new’ , f i i e=f name, access=’ sequential ’ , 

1 reel =80, form=’ formatted’ ) 
c write the occultation data to the *.oc file 

write (9, 251) i dayyr , i month 
251 format (2i5) 

wri te (9, 252) t start , tbl ock, pi tch, yaw,rol 1 , ti nc 
252 format (4d20. 12) 

write (9, 253) data 
253 format (Bf 10.7) 
close (9) 
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subroutine rol 1 i t (xsat, rsat , slat , si on, pitch, yaw, rol 1 , 
euni t, sol ax s, rl , r 2, r3, sunzen) 

this routine prepares the rotation matrices 
for the coordinate transf or mat ions required 

real *8 xsat (3) , rsat, si at, sJ on, pi tch, yaw, rol 1 , cosl , cos2, 
sinl , sin2, sunzen, phi 4, phi sun, sol ax s (3) 
real *8 r 1 (3, 3) , r2 (3, 3) , euni t (3),xunit(3),ri:(3,3), yuni t (3) 
use the following convention for the sines and cosines 
one = 1 at : two = long 

cosl = decs (si at) 
sinl = dsin(slat) 
cos2 = dcos(slon) 
sin2 = dsin(slon) 

set up lat-lon rotation matrix rl 


rl 

(1, 

r 

1) 



sinl*cos2 

rl 

(1, 

2) 

= 

sinl*sin2 

• rl 

(1, 

3) 

= 

-cosl 

rl 

(2, 

1) 

= 

-sin2 

rl 

(2, 

2) 

= 

cos2 

rl 

(2, 

3) 

= 

O.OdOO 

rl 

(3, 

1) 

= 

+cosl icos2 

rl 

(3, 

2) 

= 

+cosl*sin2 

rl 

(3, 

3) 

- 

sinl 


use now the following convention 

one = sol ar zenith angle : two = azimuth 

transform unit vector 

call rotmtx (euni t , rl , xuni t , 1 ) 

invert x and y axes 

xuni t ( 1 ) =-xuni t ( 1 ) 

xuni t (2) =-xuni t (2) 

phi 4=datan2 (xuni t(2),xunit(l)) 

compute azimuth of solar rotation matrix 

cos2=dcos (phi 4) 

si n2=dsi n (phi 4) 

cosl = dot (eunit, xsat) /rsat 

sinl = dsqrt ( 1 . dOO-cosl **2) 

sunzen=datan2 (sinl , cosl ) 


set up zenith angle 


azimuth 


rotation matrix 


r2 


r2 

r2 

r2 

r2 

r2 

r2 


(1,1) 

S 

+cos2*cosl 

(1,2) 

= 

si n2*cosl 

(1,3) 

= 

-sinl 

(2,1) 

= 

-si n2 

(2,2) 

= 

+cos2 

(2,3) 

= 

0. OdOO 
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r 2(3,1) = sinl*cos2 
r2(3,2) — sinl*sin2 
r 2(3,3) — cos 1 

compute solar axis in satelite system 

cal 1 rotmtx (sol axs, rl , xuni t , +1 ) 

x un i t ( 1 ) =->: un i t ( 1 ) 

xuni t (2) =-xuni t (2) 

call rotmtx (xuni t, r2, yuni t, +1 ) 

phi sun=datan2 (yuni t (2) , yuni t ( 1) ) 

cosl=dcos (phi sun) 

sinl=dsi n (phi sun) 

now compute the actual rotation matrix r3 
r3 (1,1) =cosl 
r3 ( 1 , 2) =sinl 
r3 (1,3) =0. OdOO 
r3 (2,1) =-si n 1 
r3 (2, 2) =cosl 
r3 (2, 3) =0. OdOO 
r 3 ( 3 , 1 ) =0. OdOO 
r3 (3, 2) =0. OdOO 
.r3 (3, 3) =1 . OdOO 
return 
end 
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SUBROUTINE ROTMTX 


ORIGINAL PASL ® 
OF POOR QUALl* / 


c 

c 

c 

c 

c 


5 


10 

15 


20 


subroutine rotmtx (xin, rot, xout, i-flag) 

this routine trans-forms column vector xin to xout, 
under the rotation matrix rot..... 
i-f i-flag less than zero the transpose o-f rot is used 

real*B xin (3) , rot (3, 3) , xout (3) 

do 5 i=l,3 

xout(i) = O.OdOO 
continue 

if( iflag .It. 0 ) goto 15 
do 10 i = 1 , 3 
do 10 j=l,3 

xout(i) = xout(i) + rot (i , j) *xin ( j) 
continue 
return 
continue 
do 20 i=l,3 
do 20 j = l , 3 

xout(i) = xout(i) + rot ( j, i ) *xin ( j) 
continue - 
return 
end 


SUBROUTINE SOLEPH 


ORIGINAL PAGE ?S 
0F Po °R QUALITY 


sub rout i ne sol eph ( ti me, dayno, :<sun, oecl ) 
c 

c this routine computes the solar ephemeris 

c -for the year 1980 using the U.S. Naval Observatory 

c chebyshev ■fit published in the Aids for Computers, 

c output is in earth-centered cartesian system 

c 


real *8 ti me, x sun (3) , crtasc (25) , cdec (25) , cdi st (25) 
real *8 

brtasc (25) , bdec (25) , bdi st (25) , :< , pi , si , rtasc, decl , distau 
real *B solrad 
integer*2 dayno 

c computed solar geocentric x,y,z 

c for dayno day of year 1980 — input 

c for time in universal time — input 

c gives :<5un in km 

data crtasc/ 

1 61 . 379B9B4d00, 1 1 . BBB27B4d00, 0. 0274567d00, 0. 0728274d00, 

2 0. O3993B0U00, 0. 1 047579d00, -0. 0307469d00, -0. 044 1 609d00, 

3 0. 0071710d00, 0. 0074422d00, -0. 0020151b00, -0. 0017746d00, 

4 O. 0010209d00, 0. 0006BlBd00, -0. 0003326d00, -0. 0001905d00, 

5 0. 00013.04d.00, 0. 00G06B6d00, -0. 0000403d00, 0. 0000069600, 

6 -0. OOOOOBSdOO, 0. 0000025d00, 0. OdOO, 0. OdOO, 0. OdOO/ 
data cdec/ 

1 -13. 4B5026d00, -2. 371 409d00, -22. 51 1744d00, 2. 734964d00, 

2 6. 654348d00, -0. 370476d00, -0. 499175d00. 0. 064997d00, 

3 0. 073062d00, -0 . 0386 1 2d00 . -0 . 036087d00, 0.01 2058d00 , 

4 0. 00832 ldOO, -0. 003071d00, -0. 001719d00, 0. 001055d00, 

5 0. 00051 OdOO, -0. 000433d00, -0. 000207d00, 0. 000130d00, 

6 -0. 00001 4d00 , 0. 000032d00, 0. OaOO, 0. OdOO, 0. OdOO/ 

data cdist/ 

1 

1 . 9B998373d00, 0. 000400B2dOC, -0. 01629251d00, -0 . 00045237d00, 

2 0. 00500746d00, 0. 00006182d00, -0. 00042170d00, 0. 00000621d00, 

3 0. 00000367d00, 0. 00000513d00, 0. 00000209d00, 0. 0000056Bd00, 

4 0. 00000356d00, 0. 00000282d00. 0. 0000052Bd00, -0. 00000141d00, 

5 


0. 00000392d 00 , -0. 00000569d00, 0. 00000007d00, -0. 00000566d00, 
6 -0. 00000470d00, 0. OdOO, 0. OdOO, 0. OdOO, 0. OdOO/ 

c 
c 


x= (dble (dayno) +ti me/24. d00+. 591 d-03) / 183. dOO-1 . 00554644Bd00 
x=2. d0*:< 

pri nt#, ’ sol eph :<= ? ,x 
do 200 i =25, 1,-1 
print*, ' i= , i 
if (i . 1 1. 23) go to 100 
brtasc ( i ) =0. OdOO 
bdec ( i ) =0. OdOO 
bdi st (i ) =0. OdOO 
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ORIGINAL 

OF POOR QUALITY 


SUBROUTINE SOLEPH 


go to 200 
100 continue 

brtasc (i ) =x*brtasc (i +1 ) -brtasc (i+2) +crtasc (i ) 
bdec(i) - x*bdec (i +1 ) -bdec (i +2) +cdec (i ) 
bdist(i) “ x*bdist (i+1 ) -bdi st (i +2) +cdist (i ) 
200 continue 

print*, ’ brtasc’ , brtasc 

print*,' bdec’ , bdec 

print*,’ bdist’, bdist 

rtasc = 7. 5d00* (brtasc ( 1 ) -brtasc (3) ) 

if (rtasc. gt. 300. dOO) rtasc = rt.ssc-360.d00 

decl = . SdOO* (bdec (1 )-baec (3) ) 

distau = . 5d00* (bdist ( 1 ) -bdist <3) ) 

sol rad = 1 . 4959B5d0B*di stau 

pi = 4. d00*datan ( 1 . dOO) 

sl= sol rad* (dcos (pi *decl / 180. dOO) ) 

xsun(l) = si *dcos (pi *rtasc/lB0. dOO) 

xsun(2) = si *dsi n (pi *rtasc/lB0. dOO) 

xsun(3) = sol rad*dsi n (pi -tdecl /1B0. dOO) 

print*,’ in soleph :;sur.= ’ , x sun 

return 

■end 


FUNCTION SEVAL 


r 

ORiGircw. r-r** r* 

OF POOR QUALITY 


real *8 -function seval (n, u, :: , y , b, c, d) 
integer n 

real *8 u, x (n),y(n),b(n),c(n),d(n) 

ccc 

ccc this subroutine evaluates the cubic spline -function 
cc 

ccc seval = y(i) + b (i > * (u-x (i ) > + c (i ) * (u-x (i ) ) **2 + 

d (i ) * (u-x (i ) ) **3 

ccc 

ccc where x(i> .It. u „lt.x(i+l), using horner's ruie 
ccc 

ccc if u .It. x (1) then i = is used 
ccc if u .gt. x (n) then i =n is used 
ccc 

ccc input, 
ccc 

ccc n = the number of data points 

ccc u = the abcissa at which the spline is to be eval . at od 
cc x,y = the arrays of data abcissas e.nd ordinates 
ccc b,c,d = arrays of spline coefficients computed b. 
cc 

ccc if u is not in the same interval as the previous .all. J 

ccc binary searcTi is performed to determine the prepo* • r " »•* • * ; - 
ccc 


int 

eger :,j> 

k 




rcu 

. *8 dx 





H -»t 

„ i/1/ 





if 

(i .ge. n 

) 

i = 

1 


i f 

( u .It. X 

(i 

) ) 

go 

to 10 

if 

(u . le. x 

(i 

+ 1 

) ) 

go to 30 


ccc 

ccc binary search 
ccc 

10 i = 1 

j = n+1 


k = 

<i+j)/2 



if 

(u .It. 

x <k> ) 

J 

= k 

if 

(u .ge. 

x ( k ) ) 

i = 

=k 

if 

(j .gt. 

i+1) 

go 

to 20 


ccc 

ccc evaluate spline 
ccc 

30 dx = u - x ( i ) 

seval = y (i ) + dx*(b(i) + dx*(c(i> ♦ 
return 
end 
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SUBROUTINE LATLON 


0f Poo fi I, 


OtJtLlf? 


subroutine 1 at Ion (xsat, t, rsat , si at , si on) 
c 

c routine computes radius, latitude and longitude 

c giver cartesian coordinates as input 

c 

real *B xsat (3) , si at, si on, rsat , rxy, t 

ruat = xsat(l)**2 + xsat(2)**2 

rxy = dsqrttrsat) 

rsat = dsqrt (rsat+xsat (3) **2) 

slat «= datan2 (xsat (3) ,rxy) 

slon = datan2 (xsat (2) , xsat ( 1 ) ) 

return 

end 
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SUBROUTINE LATLON 


ORIGINAL l» AGE 19 
OF POOR QUALITY 


subroutine 1 at 1 on (xsat, t , rsat , si at, si on) 

routine computes radius, latitude and longitude 
given cartesian coordinates as input 

real *S :<sat (3) , si at, si on, rsat , rxy, t 

rsat = xsatil)**2 + xsat(2)**2 

rxy ■ dsqrt(rsat) 

rsat «* dsqrt (rsat+xsat (3) #*2) 

slat = datan2 (xsat (3) , rxy) 

slon = datan2(::sat (2) ,xsat ( 1) ) 

return 

end 


SUBROUTINE ELL I PS 


ORIGIN /U PAGE 59 
OF POOR QUALITY 


subroutine el 1 ips (x tan, t , rtan, htan, tanl at , tan Ion, 
sundec, sunzen, pi > 

computes elliptical earth position 

implicit real *8 (a-h,o-z) 
real <4 htan 

cal 1 latl on (xtan, t, rtan, vnl at „ tanl on) 
coshr= ( (dcos (sunzen) ) -dsin (tanlat) * (dsin (sundec) ) ) / 
(dcos (tanl at) * (dcos (sundec) ) ) 
si nhr*=sqrt ( 1 . -cosi:r#*2) 
hr«*datan2 (sinhr , coshr) 
tanl on® ( ( t/24. dOO) *2*pi -pi ) +hr 

re=637B. 38B* ( 1 . 3367003e-2*dsi n (tanl at ) **2+. 70B5e-5* 
dsi n (tanl at*2) **2) 
htan=rtan-re 
tarlat=tanlat*57. 2«?r.7B 
tanl on=tanl on*57. 2957B 
return 
end 
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FUNCTION DOT 


ORfGIAiAL PAGE iS 
F POOR QUALITY 


double precision -function dot(;;,y) 
computes dot product o-f vectors x and y 
returns as dot 

real x (3) , y (3) 

dot = x(l)*y(l) + x(2)*y(2> + x(3)*y(3) 

return 

end 


APPENDIX B 


This appendix forms a brief user's manual for the programs described 
in the remainder of the report. It is intended to provide a user with all 
information necessary to be able to run the programs and obtain an inver- 

* sion of the UVSP data for a single absorbing species such as ozone. 

The first step is to get the SMM data into the computer in a usable 
format. The data was provided to Visidyne on magnetic tape. Our procedure 
was to read the files on this tape, copy them onto disk where they are much 
f simpler to handle. If these programs are to be run on the UVSP data com- 

puters then this is already accomplished, otherwise the user must get the 
data into the main disk space or directory of the target machine. Each 
UVSP experiment has been assigned a number in the form "Vxxxxx" where V 

* identifies it as a UVSP experiment and xxxxx represents the numerical order 

( of the particular experiment. (For example, the experiment number V00512 

was the very first one done for the purpose of ozone inversion.) It is 
convenient to name each file with its corresponding number. The UVSP 
"final data" files' are therefore named with the usual DEC extension as 
"Vxxxxx. FD" with "FD" standing for "final data". 

The second step is to obtain the accurate SMM ephemeris data from the 
satellite ephemeris group at Goddard Space Flight Center. This too will be 
provided on tape and the data must also be read and copies onto disk. The 
ephemeris data files are named "Vxxxxx. AO", following the convention 
described above, Note that it is necessary to be sure that the data is in 
a format appropriate for the target computer. It may be provided in one of 
several manufacturer's formats and it is best to avoid having to program 
conversions from one floating point binary format to another. 

It is useful to point out once again that the two different time 
series involved (the final data "Vxxxxx. FD" and the ephemeris "Vxxxxx. AO") 
must be accurately matched in order to obtain a valid inversion. The 
actual occultation occurs over a period of only a few seconds so it is 
clear that this can be critical. In fact several occultations could not be 
analyzed because the two series could not be matched to give reasonable 
ozone profiles at all. 

Having obtained the data as described above, the user must first run 
the program DATFND to identify which blocks of the *.FD files contain the 
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occultation. While it would certainly be possible to program the selection 
of these blocks, we chose not do so because of the relative number of 
experiments to be analyzed. It is necessary jto note which blocks do 
contain the data for each occultation experiment to be analyzed and record 
them for later use. 

The program DATPGM will prompt the user for the block numbers recorded 
from running DATFND. After this DATPGM will produce the shortened data 
file "Vxxxx.OC" and will print a set of parameters characterizing the 
experiment including such values as the wavelength and the start and end 
times of the record containing the occultation blocks. It is useful to 
keep this printout in order to know what values of the times to feed into 
the ephemeris program AOTOPT which will prompt for them. After reading 
these values, AOTOPT reads the long ephemeris file "Vxxxxx.AO" and writes 
to the shortened ephemeris file "Vxxxxx.PT", fifty satellite positions 
corresponding to a period of 500 seconds of time containing the occultation 
blocks within it. This was done for convenience in the use of the spline 
fitting and evaluation routines but a number other than fifty could just as 
well have been chosen as long as there were sufficient accuracy in the 
resulting spline fits. 

The program TANSMM uses both the "Vxxxxx.OC" and the "Vxxxxx.PT" files 
to produce a file of tangent heights corresponding to the instantaneous 
lines-of-sight for each of the data points in the two occultation blocks 
(512). The user is prompted for the experiment name and for the cross 
section for the absorbing species at the wavelength of the instrument for 
this experiment (as printed out by program DATPGM). This data is written 
out to a file called "Vxxxxx.TH". This latter file is used by the binning 
and smoothing program "BINS" along with the file "Vxxxxx.OC" to produce 
another file of a size convenient for the inversion program (REALOZ) to 
analyze. This file is called "Vxxxxx.BN". 

The reader will probably have noticed that the procedure described 
above is interactive and requires frequent responses from a user. This 
method of procedi ng was found more appropriate than a pure batch type of 
operation. Also, the process was broken up into several relatively small 
programs to simplify the use of a small computer. Most small systems do 
however provide for the definition of strings of commands which will 
eliminate the necessity of keying in these commands to run the individual 
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progams. We give two examples of such command strings, the first for the 
UNIX system and the second for the RT-11 operating system. 

UNIX SHELL COMMAND 
datpgm 
aotopt 
tansmm 
bins 
realoz 

RT-11 COMMAND STRING INTERPRETER 
R DATPGM 
R AOTOPT 
R TANSMM 
R BINS 
R REALOZ 

The UNIX commands which we shall refer to as a "shell script" are to 
be written to a file using one of the editors and the file is to be 
declared as executable after it has been created. The shell script has 

been written with the assumption that the data files reside on the same 
directory as the programs. This could be easily rearranged for another 
directory structure by anyone familiar with the UNIX shell. Note that the 
shell script does not provide for the removal of the intermediate files 
produced by the procedure. If there is a shortage of disk space, one might 
want to add a command to remove these files from the disk before finishing. 
Also one might want to provide for the redirection of the REALOZ output 
either to a printer or to a file for later printing. As the programs 
stand, the output goes to the terminal. 

The RT-11 CSI (Comnand String Interpreter) file is to be created by an 
editor and saved as a ".COM" file. (Other DEC operating systems have 
similar structures but the names and extensions may be different.) The 
procedure as written assumes that the executable files for the programs 
exist on the system device "SY :" and that the data resides on the default 
device "DK:". The output of the programs goes to the DEC line printer 
device. If the target system does not have a line printer, it will be 


necessary to assign the line printer device to the terminal or to change 
the programs to "TYPE" the output instead of "PRINT". In any case, the 
prompts currently use the "PRINT" statement and this will have to be 
changed to TYPE if one wants to use the terminal and the line printer at 
the same time. Also, the "READ" statements for terminal input will have to 
be changed to "ACCEPT" statements unless they are redirected through 
appropriate ASSIGN commands in RT-11. There are also some minor 
differences between the UNIX and DEC forms of the Fortran OPEN statement. 
It may be necessary to make some modifications to these statements if RT-11 
Fortran IV is to be used. If any of the DEC Fortran-77 compilers is used, 
however, there should not be any necessity for changing the OPEN statements 
which correspond to the ANSI standard for Fortran-77. 


APPENDIX C: LINE-OF-SIGHT GEOMETRY FOR UVSP XCULTATION EXPERIMENTS 


The inversion of the UVSP occulation results requires precise deter- 
mination of the pointing of the instrument during the occultation. Assuming, 
that the satellite position as a function of time is known to the required 
degree of accuracy, one may use the pitch, yaw and roll angles as received 
from the telemetry stream to derive the necessary pointing data. These angles 
are with respect to the instantaneous line- of- sight to the center of the solar 
disk and the solar axis. 

It is most convenient to work in the familiar Cartesian geocentric 
corrdinate system defined by the x-axis pointing at the first Point of Aries 
and the positive z-axis through the North pole. We shall call this the GCS 
system. The solar coordinates are obtained from a 22 term Chebyshev fit as 
published by the Naval Observatory. The satellite ephemeris is provided on 
tape in the GCS system as well. The unit vector of the reference line of 
sight to the center of the solar disk is given by 


(S sun - * r sat ,/l!i s U n 



U) 


where r sup and r sat are the vectors defining respectively the solar and satel- 
lite positions in the GCS system. 

A 

Any point x along the line of sight e Q may be represented in para- 
metric form the relation 


r = 


sat 


+ t e. 


( 2 ) 


where t is the parameter. It is a simple matter to determine the value of t 
which corresponds to the tangent point. Note that this t is the magnitude of 
the distance from the satellite to the tangent point. With this value of t 
one can easily compute the vector of the tangent point along the e Q line of 
sight, x tfln using Equation (2). 

We now need to determine the pointing of the sensor more precisely. 
The line of sight is specified by three angles, pitch ( 5 ), yaw ( n) and roll 
(u) with respect to the projection of the solar pole and equator on the solar 
disk. This can be handled more simply in a coordinate system with origin at 
the satellite. To transform to the new system we define a set of rotation 
matrices. First a rotation through * sat followed by a rotation through (»/2 - 

x sat> where 
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(3) 


and 


♦sat ■ ta "‘ (y sat /x sat> 


‘sat • tan '‘ ^sat /[x sat 2 * *sat 2]1/2 > 
The rotation matrix is 



sat 


sat 


0 

1 

0 



sat 

sat 

0 


sin* 

cos* 


sat 

sat 


(4) 


si "‘sat c05 *sat 

sin ‘sat slnt sat 

_cos ‘sat 

- 5iM sat 

cos *sat 

0 

+cos x sat cos» sat 

♦cos » sat sin» sat 

s,m sat 


(5) 


This is followed by a translation up the radius to the satellite by 
the distance 


sat 


sat 


+ y 


sat 


+ z 


2 1/2 


sat 


( 6 ) 


The next part of the transformation is an inversion of the transformed (x-y) 
plane which is equivalent to a rotation through * about the transformed z 
axis. The new y-axis is locally horizontal at the satellite. This is 
followed by a rotation about this latter y-axis through the zenith angle of 
the sun at the satellite. This angle is 


t = cos" (e 0 . x sat /h sat 
These two rotations are combined as 


) 


(7) 



0 

1 

0 



0 

-1 

0 



0 

-1 

0 


■sin? 

0 

cost 


) 


( 8 ) 


We now have a coordinate system whose z-axis points toward the center 
of the disk and whose y axis is really horizontal along the line of sight and 
whose x axis is parallel to the positive vertical at the tangent height. We 
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shall erect -a plane at the tangent point perpendicular to the line of sight 
and project the disk onto this plane. To define the axis on this projected 
plane we need to know the direction of the solar rotation axis on this plane. 

This direction is defined by three quantities: the angle of the 

ecliptic plane with the equator e = 23° 26', the tilt of the polar axis with 
respect to the ecliptic a = 7° 15' and the longitude of the ascending node of 
the solar equator on the ecliptic plane. This latter quantity is defined with 
a secular term to account for the earth's nutation and other minor variations 

e = 73° 40' + 50.25" t (9) 

where t is time in years after 1850. 

What we need is x» the angle that the solar axis makes with the z-axis 
of the CCS. For any given value of t in Equation (9) this will be constant. 
Clearly for one occultation experiment, t is also effectively constant and 
thus we shall ignore any changes in this angle x as well. It is a simple 
matter to derive the transformation from the CCS to the solar rotation system 
as the set of rotations 

( COSB 
0 

sinB 

Using this transformation it is simple to show that the angle x is defined by 



cosx = cose cose + sina sine sine (11) 

and the unit vector in the GCS system in the direction of the solar rotation 
axis is 


U sun 


SinB cosa 
SinB sina cose - COSB sine 
sinB sina sine + cose cose 


(12) 


Let U $ be the unit vector defining the solar rotation axis in the GCS. Then 
this vector may be transformed to the satellite system with z axis pointing at 


the center of the solar disk with the axis locally horizontal at the satel- 
lite. The tangent point along this path is readily computed and one can 
easily define a great circle plane between the satellite and the tangent point 

of the solar center line of sight. We now erect a vertical plane at the 

/ 

tangent point which, by definition, is perpendicular to the line of sight 
(z-axls). The unit vector IL may be transformed to this system by applying 

b a 

the transformations already applied to the line of sight unit vector (e Q ), 
namely Rj and R,,. Thus 

A I A T 

U s = R 2 R 1 U s 1 ( C0SY i» cos y 2 * C0Sy J (13) 

A I II I 

The angle between the projection of U s on the (x - y ) plane and the x axis 
is just 


tan* = coSYg/cosTj (14) 

This defines the orientation of the solar rotation axis on the projected 
plane. To get the desired orientation we now define another rotation matrix 

I 

to bring the x axis in line with the projected solar axis of rotation namely 

cos* 
sin* 

0 

The actual pointing of the instrument is defined by three angles, 
pitch (0, yaw (n), and roll (1). The pitch is positive toward the south 
solar pole, the yaw is positive toward the solar east and the roll is clock- 
wise from the projected solar north pole. For a given pitch (t), Yaw ( n) 
combination the center of the slit is located at 

il II 

(x ,y ) = D t (-sint, - tann) (16) 

where is the distance from the satellite to the tangent point of reference. 
For the narrow aeronomy slit (V, no. 20), we assume It is accurately repre- 
sented by a line (1" x 180"). The end points of this line are then given by 

(x"(±), y"±) = (x"±D t tan (90") slnu, y"±D t tan(90") cosu) (17) 

One can easily check whether either end point is beyond the defined edge of 
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the disk by- adding into the solar ephemeris, the expression for the apparent 
solar diameter. 

For simplicity of discussion let us calculate the center of the 
projected slit (x", y"). The other aspects are computed the same way as 
needed. The 3-vector to the center of the slit is (x", y",(x" 2 + 0 t 2 ) 1 ^ 2 . 
This then defines a unit vector e Q " for the line of sight of the measurement. 
This can now be transformed to the GCS by the series of transformations 

e = Rj T R 2 T R ax T e (18) 

and this unit vector can be used to compute the tangent point for the true 
line of sight. 

We now turn to the calculation of the latitude and longitude of the 
satellite and the tangent point. The observation is made at a given timme of 
day on a given dav of the year. The latitude of any point (x,y,z) is 

tanx = z/(x 2 + y 2 ) 1/2 (19) 

To compute the longitude of the point, one may use the relation 

costj = sinfi sinx + coso cosx cosh (20) 

where is zenith angle of the sun, 6 is the solar declination, x is the 
laJtude of the point and h is the hour angle with respect to the sun at the 
point. This can be solved for h to within a sign ambiguity. The sign is 
positive for sunset and negative for sunrise. The hour angle of the point is 

h = cos" 1 [(cos^ - sin« sinx) / cos« cosx] (21) 

The longitude is then 



1 = (UT/24)*2ir-ir + h 


( 22 ) 
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NAME 


STOP 


OPEN 

"* . AO" 


1 NPUT 

START AND 

END TIME 

READ 

EPHEMERIS 

HEADER 

PRINT 

HEADER DATA 

Ot-’EN 


FI LE/ •>* . PT" 


READ EPHEMERIS RECORD 


/Time 

WITHIN 

RANGE 


PRINT EPHEMERIS DATA 
WRITE EVERY TENTH POINT TO 
"*.PT" FILE 


DONE \ 
SO TIMES? 


CLOSE FILES 

GET NEW FILE NAME 


PROGRAM AOTOPT 












original pagl rs 

OF POOR QUALITY 


INPUT FILE NAME 
OPEN "*.OC" 

READ • .OC FOR TIMES AND ANGLES 
CLOSE *.OC 


OPEN 

* .PT 

EPHEMERI 8 FILE 

CALL 

EPHTST TO 

DO SPLINE FIT 

AND 

CLOSE 

* .PT 

FILE 


CALCULATE SOLAR AXIS VECTOR 
INITIALIZE LOOP VARIABLES 
COMPUTE TIME OFFSET AS NEEDED 


LOOP ODER DATA POINTS 


COMPUTE SOLAR EPHEMCRIS ( SOLEPH ) 

DO SPLINE FIT ( SEVAL ) 

COMPUTE SATELLITE LAT-LONG (LATLON) 


GET UT <T VECTOR TO SUN CENTER (UNIVEC) 
GET ROTATION MATRICES (ROLLIT) 

CALCULATE TANGENT HEIGHT ( TATMGT ) 

PERFORM ROTATION AND TRANSLATION ( ROTMTX) 
COMPUTE TANGENT POINT OF LINE OF SIGHT 
(TAhHGT) 

COMPUTE ALTITUDE OF TANGENV POINT 
PRINT EVERY TENTH POINT 


OPEN * . TH FILE 
WRITE TANGENT HEIGHTS 
CLOSE FILE 


u 


PROGRAM TANSMM 













ORIGINAL PAGE 

OF POOR QUALITY 


INPUT FILE NAME "»« 

OPEN "*.TH" 

READ # OF POINTS AND LAT.-LONG. 
INPUT DAWN OR DUSK 


READ TANGENT HEIGHT ARRAY 
FROM * . TH REVERSING IF 
NEEDED 

CLOSE "*.TH" 



OPEN »*.OC» 

READ DATA AND REVERSE IF 
NECESSARY 
CLOSE "*.OC" 


INPUT#OF POINTS PER BIN 
INPUT # OF POINTS TO SKIP AT START 
» I NPUT # OF POI NTS TO SKI P AT END 
OPEN "».BN» FILE 
WRITE LAT - LONG TO * . BN 


LOOP TO PERFORM REQUESTED 
BINNING AND SMOOTHING 
WRITE TO *.BN 


WHEN DONE CLOSE "».BN» FILE 


STOP 


PROGRAM BINS 























